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Non-Stationary  Subdivision  for  Inhomogeneous 
Order  Differential  Equations 


Joe  Warren  and  Henrik  Weimer 


Abstract.  This  paper  provides  a methodology  for  the  systematic  deriva- 
tion of  subdivision  schemes  that  model  solutions  to  inhomogeneous  order 
linear  differential  equations.  In  previous  work,  we  showed  that  subdivi- 
sion can  be  used  to  capture  very  efficiently  the  solutions  of  homogeneous 
order,  linear  differential  equations.  The  resulting  subdivision  masks  are 
stationary  and  can  be  precomputed,  allowing  for  very  simple  and  fast  ap- 
plication of  these  schemes.  In  this  paper,  we  show  that  this  method  can 
be  extended  to  express  solutions  of  systems  of  inhomogeneous  order,  lin- 
ear differential  equations.  Even  though  the  resulting  subdivision  masks 
may  be  non-stationary,  the  masks  can  again  be  precomputed.  Thus,  the 
resulting  subdivision  schemes  capture  very  efficiently  solutions  of  inhomo- 
geneous order,  linear  partial  differential  equations. 


§1.  Subdivision  for  the  Modeling  of  Shapes 

Subdivision  is  a popular  and  efficient  method  for  modeling  shapes.  In  par- 
ticular, subdivision  describes  a continuous  shape  p as  the  limit  of  a sequence 
Pk,  k > 0 of  discrete  shapes, 


lim  pk  = p. 

K—*00 

The  beauty  of  subdivision  lies  in  the  fact  that  these  discrete  shapes  Pk 
are  linked  by  a simple  linear  transformation  S which  is  based  on  splitting  and 
averaging, 

Pk  = Sk-iPk-i- 

Figure  1 shows  an  example  of  a subdivision  scheme.  Starting  from  the 
coarse  shape  po  on  the  left,  application  of  the  subdivision  matrix  So  yields  the 
denser  shape  p\ . As  we  continue  the  process,  the  sequence  of  discrete  shapes 
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PO  Pi  P2 


Fig.  1.  Subdivision  models  a shape  as  the  limit  of  a sequence  of  discrete  shapes. 

converges  rapidly  to  a continuous  shape  p that  follows  the  original  coarsest 
shape  po  and  whose  properties  are  determined  by  the  subdivision  matrix. 

Subdivision’s  popularity  for  the  modeling  of  curves  is  due  to  the  algo- 
rithms by  Chaikin  [9],  and  Lane  and  Riesenfeld  [6].  The  breakthrough  for  the 
modeling  of  surfaces  via  subdivision  was  marked  by  the  papers  by  Catmull 
and  Clark  [2]  and  by  Doo  and  Sabin  [3].  A popular  subdivision  scheme  for 
modeling  with  triangular  meshes  has  been  proposed  by  Loop  [7],  which  was 
also  used  for  creating  Figure  1. 

§2.  Shape  Modeling  through  Differential  Equations 

Alternatively,  shapes  can  be  characterized  as  solutions  to  partial  differential 
equations.  For  example,  any  polynomial  spline  p[x]  of  degree  m satisfies  the 
differential  equation  p(m+1)[x]  = 0,  requiring  the  ( m + l)st  derivative  of  the 
spline  to  be  zero  everywhere  except  at  a fixed  number  of  knots  [1].  Other  ex- 
amples of  shapes  based  on  partial  differential  equations  are  the  polyharmonic 
surfaces,  including  Thin  Plate  Splines,  as  well  as  many  different  classes  of  fluid 
flows. 

When  modeling  with  differential  equations,  we  determine  a continuous 
shape  p that  is  a solution  to  a set  of  partial  differential  equations 

D p = 6,  (1) 

where  D denotes  a continuous  differential  operator  and  6 encodes  the  bound- 
ary conditions  for  the  problem.  For  the  example  of  natural  cubic  splines,  we 
have  D = and  6 = 0 almost  everywhere.  If  all  differential  operators  in 
D are  of  the  same,  fixed  order,  we  call  the  differential  equation  homogeneous 
order.  Otherwise,  the  equation  is  called  inhomogeneous  order. 

To  handle  such  problems  in  a computational  environment,  one  commonly 
discretizes  the  continuous  problem.  To  this  end,  a domain  grid  is  chosen 
and  all  entities  of  the  continuous  partial  differential  equation  (1)  are  discretized 
over  this  domain  grid.  The  result  is  a system  of  linear  equations 


DkPk  = bit. 


(2) 
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where  pk  denotes  an  approximation  of  the  continuous  solution  p over  the  grid 
Tk , bk  denotes  a discretization  of  the  boundary  conditions,  and  Dk  is  a discrete 
approximation  of  the  continuous  differential  operators  D on  the  domain  grid 
Tk. 

Relying  on  the  theory  for  finite  elements  or  finite  differences  [11],  the 
discrete  solutions  pk  can  be  formally  guaranteed  to  converge  to  the  continuous 
solution  p of  the  original  continuous  problem  (1)  if  the  discretizations  Tk  are 
chosen  carefully  and  the  discrete  representations  Dk  and  bk  are  well  chosen. 

At  this  point,  the  problem  of  finding  the  continuous  solution  p of  the 
system  of  continuous  partial  differential  equations  (1)  has  been  reduced  to  the 
problem  of  solving  denser  an  denser  systems  of  linear  equations  (2). 

The  links  between  mesh  modeling  and  differential  equations  were  previ- 
ously investigated  by  Mallet  [8],  Taubin  [12],  and  Kobbelt  [5].  The  method 
presented  here  is  new  because  subdivision  schemes  that  model  solutions  of  in- 
homogeneous order  differential  equations  are  precomputed  entirely,  enabling 
very  efficient  modeling  of  shapes  guided  by  inhomogeneous  order  differential 
equations.  In  particular,  the  actual  application  of  the  subdivision  schemes 
does  not  require  any  computational  solving  whatsoever. 

§3.  Subdivision  for  Homogeneous  Order  Differential  Equations 

In  our  previous  work  [13,14]  we  characterized  subdivision  schemes  for  the 
solutions  of  homogeneous  order  linear  partial  differential  equations.  In  this 
framework,  the  subdivision  matrix  Sk~i  is  determined  as  the  solution  to  the 
system  of  linear  equations 

DkSk-i  = 2dUk-\Dk-\,  (3) 

where  d is  the  dimension  of  the  domain.  Recall  that  the  differencing  operator 
Dk  is  the  discrete  approximation  of  the  continuous  differential  operator  D of 
the  original,  continuous  problem  (1)  on  the  level  k grid  Tk.  Further,  Uk-i 
denotes  a very  simple  linear  transformation,  called  replication  or  upsampling, 
that  carries  coefficients  over  the  grid  Tk-i  into  coefficients  over  the  next  denser 
grid  Tfc.  The  action  of  Uk-i  is  very  simple:  Coefficients  centered  over  knots  in 
Tk- 1 are  replicated  over  the  same  knots  in  the  denser  grid  Tk  while  coefficients 
centered  over  the  remaining  knots  Tk  — Tk-\  are  set  to  zero.  Thus,  Uk-i  is  a 
matrix  whose  rows  are  either  zero  or  a standard  unit  vector,  and  Uk- i can  be 
constructed  easily  and  efficiently. 

We  visualize  the  meaning  of  equation  (3)  in  Figure  2:  The  subdivision 
matrix  is  determined  so  that  a certain  commutativity  relationship  holds  be- 
tween subdivision,  upsampling  and  differencing.  Differencing  coefficients  on 
the  coarse  grid  and  upsampling  those  differences  to  the  finer  grid  by  inserting 
zero  for  all  new  grid  points  (Uk-iDk-i,  the  right  hand  side  of  equation  (3)) 
should  yield  the  same  result  as  subdividing  the  coefficients  using  the  subdi- 
vision scheme  and  then  differencing  on  the  finer  grid  ( DkSk-i , the  left  hand 
side  of  equation  (3)). 
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Fig.  2.  The  subdivision  scheme  is  determined  such  that  this  commutativity  re- 
lationship holds. 


The  subdivision  matrices  Sk  are  the  only  unknowns  of  equation  (3),  and 
we  can  use  linear  algebra  to  systematically  solve  for  these  subdivision  matri- 
ces. In  our  previous  work  [13]  we  showed  that  solutions  produced  by  these 
subdivision  matrices  are  related  to  solutions  produced  by  an  interpolating  fi- 
nite element  solver  using  a simple,  fixed  change  of  basis.  Therefore,  if  the  finite 
element  solver  converges,  then  the  subdivision  solution  is  also  well  defined. 

As  a side  note,  in  our  previous  work  [13]  we  establish  that  the  right-hand 
side  of  the  system,  as  solved  by  the  subdivision  scheme  from  relation  (3),  is 
Dkpk  = Uk-iUk-2  ■ ■ ■ UoDqPo  where  pk  = Sk-iPk-i  and  p0  is  a user-given 
set  of  initial  control  coefficients.  In  other  words,  the  subdivision  scheme  leads 
to  a specific  combination  of  the  integer  shifts  of  the  Green  function  of  the 
differential  operator. 

Further,  computation  of  the  subdivision  matrix  Sk-i  based  on  relation 
(3)  requires  the  inversion  of  the  differencing  operator  D*. . Consequently,  the 
computational  work  required  for  finding  the  subdivision  matrix  is  at  least 
the  same  as  inverting  the  finite  difference  system.  However,  later  we  will  see 
that  the  subdivision  matrices  can  be  precomputed.  Thus,  in  contrast  to  a 
conventional  finite  difference  solver,  new  shapes  can  be  generated  extremely 
efficiently. 

As  an  example,  we  briefly  derive  subdivision  schemes  for  piecewise  poly- 
nomial splines.  Recall  from  deBoor  [1]  that  the  piecewise  polynomial  spline 
p[x]  of  degree  m satisfies  the  differential  equation  £)[x]m+1p[x]  = 0 where 
D[x\  denotes  the  first  derivative  in  the  variable  x. 

We  employ  generating  functions  [4]  for  concise  and  convenient  encoding 
of  discrete  coefficient  sequences.  To  this  end,  we  choose  the  domain  grids  for 
our  analysis  as  the  dilates  of  the  integer  grid  TL.  A generating  function 
Pk  [a:]  is  a power  series  that  associates  the  ith  coefficient  of  the  discrete  shape 
Pk  as  the  coefficient  of  x'.  For  example,  the  coefficient  sequence  {1,2,  3, 4,  5} 
is  represented  by  1 + 2a;  + 3a;2  + 4x3  + 5a;4. 

Recall  the  definition  of  the  first  derivative  operator, 


D[x]p[x]  = lim 
t-o 


p[s]  - p[x  + 1\ 


t 
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Substituting  t = ^ yields 


D[x\p\x ] = lim 

k—*oo 


p[x]  - p [at  + £] 
1 

2* 


Thus,  for  x £ -^TL,  the  approximation  of  the  first  derivative  is  given  by  the 
difference  between  two  adjacent  discretizations,  normalized  by  the  grid  spac- 
ing. In  terms  of  generating  functions,  this  differencing  operation  is  represented 
by  the  Laurent  polynomial 

Higher  order  derivatives  and  differences  are  obtained  by  repeated  application 
of  the  respective  continuous  or  discrete  operator. 

In  terms  of  generating  functions,  the  action  of  the  upsampling  operator 
Uk  can  be  captured  very  concisely:  The  expression  p [x2]  represents  the  up- 
sampled  coefficient  sequence  of  p[x ] as  a generating  function.  Thus,  in  our 
example  of  polynomial  splines,  the  generating  function  s*  [a:]  for  the  subdivi- 
sion scheme  satisfies 

Dk[x]m+lsk[x]  = 2Dfc_1  [x2]m+1, 


which  can  be  simplified  to 


= 2 


(Dk-i  [z2]\m+1 

V Dk[x\  ) 


Fortunately, 

Dk-i[x2}  _ 1 1 + x 
Dk[x]  2 x1/2 

i.e.  the  generating  functions  for  the  differencing  operations  on  the  level  k — 1 
and  level  k grids  divide  out  yielding  a simple  expression  independent  of  k.  As 
a result,  the  subdivision  mask  for  the  degree  m polynomial  splines  are  exactly 
the  coefficients  of 

■M&r 

Remarkably,  these  are  precisely  the  known  subdivision  schemes  for  piece- 
wise  linear  functions  (m  = 1),  Chaikin’s  algorithm  [9]  (m  = 2)  and  the 
Lane/Riesenfeld  algorithm  [6]  (m  = 3). 

Previously  we  applied  this  strategy  to  derive  subdivision  schemes  model- 
ing solutions  of  homogeneous  order  linear  differential  equations  yielding  local, 
stationary  subdivision  masks  [13,14],  In  this  paper,  we  show  that  largely  the 
same  strategy  can  be  used  to  determine  subdivision  schemes  for  inhomoge- 
neous order  linear  partial  differential  equations.  As  we  will  see,  in  this  case 
the  actual  subdivision  masks  may  depend  on  the  particular  level  of  subdivi- 
sion, i.e.  are  non-stationary.  However,  the  masks  can  still  be  precomputed  as 
a closed  form  algebraic  expression  in  the  level  of  subdivision,  which  can  then 
be  evaluated  very  efficiently  during  the  actual  application  of  the  scheme. 
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§4.  Subdivision  for  Inhomogeneous  Order  Differential  Equations 

In  this  section  we  extend  our  systematic  construction  of  subdivision  schemes  to 
handle  inhomogeneous  order  linear  partial  differential  equations.  We  consider 
the  simple  yet  interesting  problem  of  splines  in  tension  [10].  The  continuous 
spline  in  tension  p[x ] for  tension  parameter  a satisfies  the  differential  equation 

(D[z]4  - a2  D[x ]2)  p[x]  = 0,  (4) 

where  D[x]  again  represents  the  continuous  first  derivative  operator  with  re- 
spect to  the  variable  x.  Note  that  equation  (4)  incorporates  both  second  and 
fourth  derivatives  of  p[sr],  i.e.  is  inhomogeneous  order. 

Following  the  same  strategy  as  in  the  derivations  for  polynomial  splines, 
we  use  generating  functions  to  encode  the  discrete  approximation  pk  of  the 
spline  in  tension  on  grid  Tk  as  well  as  for  the  representation  of  the  differencing 
operation  Dk[x ] = Next,  we  apply  equation  (3)  to  characterize  the 

subdivision  scheme  s/tjz]  as  the  solution  to 

(•DfcN4  - a2Dk[x]2^j  sfc_i[x]  = 2 [x2]4  - a2Dk_i  [x2]2)  , (5) 


which  can  be  simplified  to 


«fc-i[x] 


2 [pk- i[x2]4  - a2Dk- i[x2]2) 
dT[x]4  - aWk [x]2 


(6) 


However,  at  this  point  we  note  that  there  is  no  simple  closed-form  ex- 
pression for  the  quotient  ' (unless  a = 0).  In  other  words, 

there  is  no  finitely-supported  subdivision  scheme  Sjt_i[x]  for  splines  in  ten- 
sion. Moreover,  the  coefficients  of  the  Laurent  series  expansion  of  the  quotient 
s;t_i[x]  depend  on  the  level  of  subdivision  k , i.e.  the  subdivision  scheme  has 
to  be  non- stationary. 

Fortunately,  due  to  the  structure  of  equation  (3)  the  coefficients  of  this 
expansion  decrease  very  rapidly  away  from  the  origin.  Thus,  we  can  approxi- 
mate the  infinite  Laurent  expansion  of  the  subdivision  mask  well  by  a locally 
supported  scheme.  To  this  end,  we  construct  the  generating  function  sjt-i[x] 
of  desired  support  symbolically  with  the  actual  coefficients  s'k_1  as  unknowns, 


n 

Sk-l[x]  = Sk-lx' 

i=  — n 


for  a user-defined  support  n.  We  then  construct  a generating  function  for  the 
residual  of  equation  (5), 


rk[x\  = ^Dfc[x]4  - a2Dk[x]2^  s*_i[x]  - 2 [x2]4  - a2Dk- 1 [x2]2) 

=Y,r*xi- 


(7) 
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Using  linear  algebra,  we  can  now  solve  for  the  unknowns  of  (7)  sym- 
bolically by  minimizing  the  least  squares  residual  of  the  coefficients  r jj. . The 
motivation  behind  our  strategy  is  to  construct  a best  solution  for  the  charac- 
teristic equation  (3)  of  given  support.  The  results  of  this  process  are  actual, 
symbolic  coefficients  for  the  local  subdivision  scheme  Sfc_i[x],  depending  on 
the  tension  parameter  a as  well  as  on  the  level  k.  As  an  example,  the  ap- 
proximation to  (6)  with  the  same  support  as  the  Lane-Riesenfeld  algorithm 
(n  = 2)  has 

22  * (693  21+1° k +891  41+4 k a2+3525 64fc  a4+399  41+2 k a6+333  4k  <*®+26  a10) 

8 (693  21+'2  k+891 4l+5  k a2+3861  256k  a4+273  23+6  k a6+675  16fc  a8+27  41+fc  a10 +7  a12) 

as  the  coefficient  for  x±2, 

(41+t+a2)  (693  1024fc+171 4k  aa+2  a2  (891  256k+219  41+3  k a2+399  16fc  a4+7a8)) 

4(693  21+12k+89141+5k  <*2+3861256k  a4  + 273  23+6fc  a6+675  16*  a8+27  41+fc  a10+7  a12) 


as  the  coefficient  for  x*1,  and  finally 

2079 41+6 *+6039 41+5 k q2  + 1755  161+2t  q4+826521+6t  q6+5235  16*  q8+21341+k  q10+56q12 
8(693  2i+12k+89141+5t  a2+3861256*:  a4+273  23+6t  a6+67516t  a8+2741  + fc  a10+7a12) 

as  the  coefficient  associated  with  x°.  Note  that  for  a = 0 these  coefficients 
exactly  reduce  to  the  subdivision  scheme  for  natural  cubic  splines  based  on 
the  Lane-Riesenfeld  algorithm. 

During  an  actual  application  of  the  subdivision  scheme,  the  user-defined 
tension  parameter  a and  the  current  level  of  subdivision  k are  substituted 
into  the  symbolic  solution  Sfc_i[x],  yielding  a simple  generating  function  in 
only  the  variable  x.  The  coefficients  of  this  generating  function  encode  the 
subdivision  masks  for  the  spline  in  tension  for  the  given  tension  parameter  a 
at  the  current  level  k.  Again,  application  of  this  subdivision  scheme  is  very 
efficient.  For  example,  given  a = 0,  the  above  expression  simplifies  to  the 
generating  function  for  natural  cubic  spline  subdivision,  independent  of  k. 


k = 1 : 0.11063  0.55261 

k = 2 : 0.12345  0.52441 

k = 3 : 0.12489  0.50733 

k = 4 : 0.12499  0.50192 

k = 5 : 0.125  0.50049 

k = 6 : 0.125  0.50012 

k = 7 : 0.125  0.50003 

k = 8 : 0.125  0.50001 

k = 9 : 0.125  0.5 

k = 10 : 0.125  0.5 

k = 11  : 0.125  0.5 
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Fig.  3.  Subdivision  masks  for  a = 1,  k = 0. 
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Fig.  4.  Splines  in  tension  for  varying  a. 

Figure  3 shows  the  actual  coefficients  of  a locally  supported  generating 
function  (n  = 2)  for  a = 1 and  k = 1, . . . , 11.  Coefficients  were  rounded  to 
five  significant  digits.  Note  that  the  coefficient  sequence  rapidly  converges  to 
the  subdivision  scheme  for  natural  cubic  splines.  Indeed,  after  a few  rounds 
of  subdivision,  a spline  in  tension  behaves  like  a natural  cubic  spline  over  a 
denser  initial  grid  with  its  initial  control  coefficients  determined  by  the  first 
few  rounds  of  subdivision. 

Figure  4 depicts  application  of  four  rounds  of  the  local  subdivision  scheme 
(support  n — 4)  for  a ranging  from  0 to  5.  The  initial  control  polygon  is  shown 
as  a thin  line,  while  the  subdivided  curve  is  shown  in  solid.  Note  that  as  a 
is  increased,  the  curve  follows  the  control  polygon  more  closely.  In  the  limit, 
a — > oo,  the  curve  is  actually  the  piecewise  linear  interpolant  of  the  initial 
control  points. 

Figure  5 shows  the  least  squares  residuals  ^T(r^)2  of  approximations  of 
different  sizes  for  a = 1 and  k = 0 (the  residual  is  largest  for  k = 0)  on  a 
logarithmic  scale.  Note  that  for  the  approximation  of  size  n = 4 the  residual 
is  already  very  small. 

At  a higher  level,  we  follow  these  steps  in  the  derivation  of  non-stationary 
subdivision  schemes  for  inhomogeneous  order  linear  partial  differential  equa- 
tions: 

Starting  from  the  continuous,  inhomogeneous  order,  linear  partial  dif- 
ferential equations  we  discretize  the  continuous  differential  operators  to  yield 
appropriate  differencing  operators  over  the  respective  subdivision  grids  T*,  . We 
then  characterize  the  subdivision  scheme  Sk-i  as  the  only  unknown  of  equation 
(3)  using  these  differencing  operators  as  well  as  simple  replication/upsampling. 
Next,  we  construct  a representation  of  the  subdivision  scheme  Sk- i in  terms 
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Fig.  5.  Residuals  of  local  approximations  with  varying  support. 

of  unknowns  and  symbolically  build  a residual  expression  representing  the  dif- 
ference between  left-hand  side  and  right-hand  side  of  equation  (3).  Finally, 
we  use  linear  algebra  to  solve  symbolically  for  the  unknowns  of  the  subdivi- 
sion scheme  Sk-i,  which  may  depend  on  the  level  of  subdivision  and  possible 
parameters  to  the  original  partial  differential  equations.  As  a result,  applica- 
tion of  the  subdivision  scheme  only  involves  instantiation  of  these  constants, 
yielding  a locally  supported,  approximating  subdivision  scheme  for  solutions 
of  the  original  inhomogeneous  order  partial  differential  equations. 

§5.  Summary  and  Conclusion 

In  this  paper,  we  showed  that  subdivision  can  be  used  to  model  solutions 
of  inhomogeneous  order  differential  equations.  Using  the  characterization  of 
the  subdivision  scheme  based  on  the  commutativity  relationship  (3),  we  can 
systematically  solve  for  these  schemes.  Even  though  the  exact  subdivision 
schemes  may  be  globally  supported,  locally  supported  schemes  approximate 
the  solution  well  enough  for  practical  purposes.  Non-stationary  schemes  can 
be  handled  using  the  same  methodology  by  allowing  the  locally  supported  sub- 
division masks  to  change  between  levels.  Because  these  subdivision  schemes 
can  be  precomputed,  the  modeling  of  solutions  of  inhomogeneous  order  linear 
partial  differential  equations  can  be  handled  very  efficiently. 

The  proposed  method  for  modeling  solutions  to  inhomogeneous  order  lin- 
ear differential  equations  is  quite  general  and  promises  to  be  useful  in  a variety 
of  applications.  First  of  all,  approximations  based  on  local  subdivision  schemes 
are  often  sufficient  for  modeling  applications.  Indeed,  the  approximate  solu- 
tions are  qualitatively  indistinguishable  from  the  exact  solution.  Second,  if  the 
accuracy  of  the  subdivision  solution  is  not  satisfactory,  the  subdivision  scheme 
can  be  used  to  produce  very  good  initial  estimates  for  more  traditional  solu- 
tion methods.  Third,  the  results  of  traditional  solution  methods  often  need  to 
be  refined  locally  for  visualization  and  analysis.  A local  subdivision  scheme 
can  be  used  to  refine  solutions  to  any  desired  accuracy  and  provide  better 
accuracy  than  traditional  polynomial  fits. 
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